/* 


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-g96.c:

   Automatically generated code implementing G96 functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional g96Functional;" to 'functionals.h'
    2. add "&g96Functional," to 'functionals.c'
    3. add "fun-g96.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
b:   1.0/137.0;
xa:  grada/(rhoa^(4/3));
xb:  gradb/(rhob^(4/3));

Exa: -rhoa^(4/3)*b*xa^(3/2);
Exb: -rhob^(4/3)*b*xb^(3/2);

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer g96_isgga(void) { return 1; } /* FIXME: detect! */
static integer g96_read(const char *conf_line);
static real g96_energy(const FunDensProp* dp);
static void g96_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void g96_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void g96_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void g96_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional G96xFunctional = {
  "G96x",       /* name */
  g96_isgga,   /* gga-corrected */
   1,
  g96_read,
  NULL,
  g96_energy,
  g96_first,
  g96_second,
  g96_third,
  g96_fourth
};
 
/* IMPLEMENTATION PART */
static integer
g96_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
g96_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(rhob,1.333333333333333);

   /* code */
    res = -0.007299270072993*t2*pow(sqrt(gradb/t2),3.0)-0.007299270072993*
        t1*pow(sqrt(grada/t1),3.0);

    return res;
}

static void
g96_first_helper(real rhoa, real grada, real *res)
{    real t1;

    t1 = sqrt(grada/pow(rhoa,1.333333333333333));

   /* code */
    res[0] = 0.014598540145985*t1*grada/rhoa-0.009732360097324*
        pow(t1,3.0)*pow(rhoa,0.333333333333333);
    res[1] = -0.010948905109489*sqrt(grada)/pow(rhoa,0.666666666666667);
}

static void
g96_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    g96_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        g96_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
g96_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4;

    t1 = sqrt(grada/pow(rhoa,1.333333333333333));
    t2 = pow(t1,3.0);
    t3 = sqrt(grada);
    t4 = 1/pow(rhoa,0.666666666666667);

   /* code */
    res[0] = 0.014598540145985*t1*grada/rhoa-0.009732360097324*
        t2*pow(rhoa,0.333333333333333);
    res[1] = -0.010948905109489*t3*t4;
    res[2] = -0.004866180048662*pow(t3,3.0)/pow(rhoa,2.666666666666667)-
        0.003244120032441*t2*t4;
    res[3] = 0.007299270072993*t3/pow(rhoa,1.666666666666667);
    res[4] = -0.005474452554745*t4/t3;

}

static void
g96_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    g96_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        g96_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
g96_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8;

    t1 = sqrt(grada/pow(rhoa,1.333333333333333));
    t2 = pow(t1,3.0);
    t3 = sqrt(grada);
    t4 = 1/pow(rhoa,0.666666666666667);
    t5 = pow(t3,3.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = 1/pow(rhoa,1.666666666666667);
    t8 = 1/t3;

   /* code */
    res[0] = 0.014598540145985*t1*grada/rhoa-0.009732360097324*
        t2*pow(rhoa,0.333333333333333);
    res[1] = -0.010948905109489*t3*t4;
    res[2] = -0.003244120032441*t2*t4-0.004866180048662*t5*
        t6;
    res[3] = 0.007299270072993*t3*t7;
    res[4] = -0.005474452554745*t8*t4;
    res[5] = 0.019464720194647*t5/pow(rhoa,3.666666666666667)+
        0.002162746688294*t2*t7;
    res[6] = -0.012165450121655*t3*t6;
    res[7] = 0.003649635036496*t8*t7;
    res[8] = 0.002737226277372*t4/t5;

}

static void
g96_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    g96_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        g96_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
g96_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;

    t1 = sqrt(grada/pow(rhoa,1.333333333333333));
    t2 = pow(t1,3.0);
    t3 = sqrt(grada);
    t4 = 1/pow(rhoa,0.666666666666667);
    t5 = pow(t3,3.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = 1/pow(rhoa,1.666666666666667);
    t8 = 1/t3;
    t9 = 1/pow(rhoa,3.666666666666667);
    t10 = 1/t5;

   /* code */
    res[0] = 0.014598540145985*t1*grada/rhoa-0.009732360097324*
        t2*pow(rhoa,0.333333333333333);
    res[1] = -0.010948905109489*t3*t4;
    res[2] = -0.003244120032441*t2*t4-0.004866180048662*t5*
        t6;
    res[3] = 0.007299270072993*t3*t7;
    res[4] = -0.005474452554745*t8*t4;
    res[5] = 0.002162746688294*t2*t7+0.019464720194647*t5*
        t9;
    res[6] = -0.012165450121655*t3*t6;
    res[7] = 0.003649635036496*t8*t7;
    res[8] = 0.002737226277372*t10*t4;
    res[9] = -0.075696134090295*t5/pow(rhoa,4.666666666666667)-
        0.003604577813824*t2*t6;
    res[10] = 0.032441200324412*t3*t9;
    res[11] = -0.006082725060827*t8*t6;
    res[12] = -0.001824817518248*t10*t7;
    res[13] = -0.004105839416058*t4/pow(t3,5.0);

}

static void
g96_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    g96_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        g96_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
